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The growth in computer hardware performance and capacity has enabled 
large scale computations of complex physical models. 
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These calculations raise several questions: 
o How accurate is the simulation? 
o Can predictions be trusted? 

o Can differences between computation and experiment be rigorously 
reconciled? 
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o The main topic of discussion is error representation and error control of 
functional outputs via dual problems (Erickson et. al., 1995), (Becker 
and Rannacher, 1997). 

o Particular attention is given to the time-dependent calculation of the 
compressible Navier-Stokes flow. Specifically, we examine the 
backwards-in-time dual problem and issues associated with 
o the deterioration (blowup) of dual problems with increasing Reynolds 
number, 

o the loss of sharpness in error bounds over long time integrations. 

o In the remainder of the presentation, we briefly examine a novel 
uncertaintly quantification technique proposed by Estep and Neckels 
(2006) for the quantification of uncertain functional outputs given 
aleatoric (statistical) random variable inputs. 

o Surprisingly, the dual problems in the Estep and Neckels technique are 
identical to those arising in a posteriori error estimation (!!) but now the 
dual problem is used to construct a piecewise linear approximation of 
the random variable response surface. 

< □ M g M I M I ► 1 -0 0,0 


Motivating Computational Challenge #1 : Cylinder Flow 



Space- 

Time 

DG 

Tim Barth 

introduction 

Cylinder 

Flow 

Nonlinear 

Conserva- 

tion 

Laws 

Space-time 

DG 

Error Rep- 
resentation 

Scalar 

transport 

Example 

Dual 

Problems 

Navier- 

Stokes 

Formulation 

Periodic 

Cylinder 

Flow 

Concluding 

Remarks 


Cylinder flow at Mach = 0.10, logarithm of |vorticity| contours 



Re=1000 



Re=3900 



Re=10000 Re=50000 


o Quartic space-time elements . | '~ /" — 

o 25K element mesh 

o Viscous walls only imposed on cylinder surface 
o Reynolds number based on cylinder diameter 



Question: How is the ability to estimate and control numerical error effected 
by increasing Reynolds number? 
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Conservation law system in R dx1 

u,t + divf = 0, u,f/eR” / = 1,...,d 


Convex entropy extension 

U,t + div F < 0, U, Fj £ R 


□ ► < 
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Space-Time Discontinuous Galerkin Formulation 
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Piecewise polynomial approximation space: 

V h = {v„ | V/,| Kxi n £ ( 'Vk(K x / n )) } 

Find \ih £ V h such that for all w h e V h 

W-1 

S(v^,w h ) D G = Yh B n (Vh:Wh)DG = 0 , 

n— o 

e n (v,w) DG = f Y' f -(u(v) • w it + f'(v) • w, x ,)£*rcff 

Jln i<£r J K 

+ / y 2 1 w(x_) ■ h(v(x_), v(x + ); n) dsdt 

Jl" K£T J 8K 

+ J (w(f" +1 ) ■ u(v(f” +1 )) - w(f”) ■ u(v(f"))^ dx 

o Proposed by Reed and Hill (1973), LeSaint and Raviart (1974) and further 
developed for conservation laws by Cockburn and Shu (1990) 

o u the conservation variables, v the symmetrization variables 
o ha numerical flux function, h(v_, v + ; n) = -h(v + , v_; -n), h(v. v; n) = f(v) • n 

< g ► <!► <!► 1 -o o 


The Discontinuous in Time Approximation Space 



Space- 

Time 

DG 

Tim Barth 

Introduction 

Cylinder 

Flow 

Nonlinear 

Conserva- 

tion 

Laws 

Space-time 

DG 

Error Rep- 
resentation 

Scalar 

transport 

Example 

Dual 

Problems 

Navier- 

Stokes 

Formulation 

Periodic 

Cylinder 

Flow 

Concluding 

Remarks 


o Natural setting for the discontinuous Galerkin (DG) method for 
hyperbolic problems 

o Utilized in the space continuous Galerkin least-squares method 
(Hughes and Shakib, 1988) 

o Often used in the discretization of parabolic problems (Douglas and 
Dupont, 1976) 

o Requires solving the implicit slab equations. 



Discontinuous timeslab 
intervals 


Space-time prism element 
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Nonlinear Stability of Space-Time DG Formulations 
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Theorem E: Global space-time entropy inequality (Cauchy IVP): 

[ U(u*(t°_))dx< [ U(u(v h (x,t N )))dx< f U(u(v h (x,t°_)))dx 
Jo. Jo. Jn 

u ‘ ( '-> = ^mL u[Mx ’ f - ))dx 

whenever the numerical flux satisfies the system extension of Osher’s 
famous “E-flux” condition 


[v]L(h(v_,v +: n)-f(v(0)).n)<O , V© e [0,1] , v(f?) = v_ + 0[v]: 


Navier- 

Stokes 

Formulation 


o Several flux functions satisfy this technical condition when recast in 
entropy variables, e.g. Lax-Friedrichs, HLLE, Roe with modifications, 
etc. 
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and Theorem E is satisfied for the Cauchy IVP, then following /_ 2 stability 
result is readily obtained 


< (c ^ 1 C 0 ) 1/2 ||u(v„(-, f° )) - u’(?° )||, 
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Space-Time Error Control 
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Given a system of PDEs with exact solution u e R m in a domain O, the overall 
objective is to construct a locally adapted discretization with numerical solution uu that 
achieves 

o Norm control [Babuska and Miller, 1984] 

llu — u/,11 < tolerance 

o Functional output control [Erickson et. at. (1995), Becker and Rannacher, 1997] 
|J(u) - J(u/,)| < tolerance , J(u) : R" 1 1 — > R 

Example functional outputs: 

o Time-averaged lift force, drag force, pitching moments 
o Average flux rates through specified surfaces 
o Weighted-average functionals of the form 

J v (u) = [ f 4r(x, t) • N(u)dx dt 
Jt 0 Jn 

for some user-specified weighting rlf(x, f) and nonlinear function N(u ) 

<!► <!► 1 -0 0,0 


Error Representation: Linear Case 
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Assume S(-, ■) bilinear and J(-) linear. 

Primal Numerical Problem: Find u^, e V® such that 

B(u h , w) = F(w) VweV B . 

Auxiliary Dual Problem: Find <t> e V B such that 

S(w. <t>) = J(w) VweV B . 

J(u) - J(Uh) = J ( u - Uf,) (linearity of J) 

= B(u - U/,, <t>) (dual problem) 

= B(u - Uh, <t> - ir h <i>) (Galerkin orthogonality) 

= B(u, O - n h <t>) - B(u h , <t> - 717,0) (linearity of B) 

= F(<t> - n h <t>) - B(u h , 0 - 77,0) (primal problem) 

Final error representation formula: 

J(u) - J{ u h ) = F(c t> - 717,0) - B( u„, <t> - n h <t>) 

< □ ► < g ► < 1 ► <!► 1 -0 0,0 
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Various techniques in use for estimating 0 - 717 , 0 : 

o Higher order solves [Becker and Rannacher, 1998],[B. and Larson, 
1999], [Suli and Houston, 2002], [Houston and Hartman, 2002] 

o Patch postprocessing techniques [Cockburn, Luskin, Shu, and S uli, 
2003] 

o Extrapolation from coarse grids 



— (u(6)) dd ■ {u B - u A ) = L u ■ (u B - U A ) 
Jo du 


B(u,w) = B(Uh, w) + (L u ■ (u - u h ),w) 

= B(u h , w) + £?(u - u h: w) W e V B 


1 'O 
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Semilinear form B{-, ■) and nonlinear J( ). 

Primal numerical problem: Find u h e such that 

B(u h , w) = F(w) VweV B . 

Linearized auxiliary dual problem: Find <t> e V B such that 
£?(w, <t>) = J(w) VweV B . 

J( u) - J(u h ) = J(u - u ft ) (mean value J) 

= H(u-u h ,<t>) (dual problem) 

= B( u - u h , <t> - 777,0) (Galerkin orthogonality) 

= B( u, O - 7 77,0) - B(Uh , 0 - -ny, 0 ) (mean value B) 

= F (0 - 777,0) - B(Uh, 0 - 717,0), (primal problem) 

Final error representation formula: 

J ( U) - J ( U h ) = F (0 - 777 , 0 ) - B{ U h , 0 - 77 ^, 0 ) 

< □ m9 m | ► < | ► 1 -0 0,0 



BdG ( v hi w ) ~ Fdq(* ~ — X X ®DG ,Q"( v h- ^ ~ 7r h < t > ) ^ F DG,0"(‘I > ~ 7 '7i < t > ) 

n— 0 Q" 


X] X S DG ,Q n ( V/7 , <t> — 7r/,<t>) — F DG,O n ( < * > ~ 7r /7 < l > ) 
A7=0 O n 


|J(u) - J(u h )| < X X l^o.Q"^. ^ - % < 1 , )| 

n = 0 Q n ' 

refinement indicator , i 


Fixed fraction mesh adaptation: 

o Refine a fixed fraction of element indicators, rj Q n, that are too large and coarsen 
a fixed fraction of element indications that are too small. 


Example: A Scalar Time-Dependent PDE 
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Circular transport, A = (y, -x), of bump data 


Uf + A- Vu = 0 , xe[-1,1] 2 

u(x, 0) = Ur(1 /10; x — x 0 ) , x 0 = (7/10,0,0) 


Primal numerical problem Convergence, ||u — z_ 2 (n x [o, 7 "]) 


□ i >00.0 





Example: A Scalar Time-Dependent PDE 


Tim Barth 
Introduction 


A functional is chosen that averages the solution data in the space-time ball of radius 
1/10 located at x c = (1 /2, 1 /2, 1 .05) in space-time 

C 1 .15 

i dxdt 

Jo Ja 


J( u) = [ [ r|r(1 /1 0; x — x c ) U ( 

Jo Ja 
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J(u) - J( U h ) = F DG,Q"(«> - 7T/7<t>) - 6D0,0"(Vft, <t> - n h < t>) 

n=N--\ K 
0 

\J(U) - J(u h )\ < -Kh®) - B DGtQ n(\l h ,<S> --rr h <!>)\ 

n=N - 1 K 




Error egimat§ buildup 
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An Application of Error Estimation and Adaptive Error Control 
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Example: Euler flow past multi-element airfoil geometry. M = .1 , 5° AOA. 
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lift coefficient 
(error representation) 

lift coefficient 
(error control) 

refinement 

level 

# elements 

equivalent uniform 
refinement 
# elements 

5.156 ± .147 

5.156 ± .346 

0 

5000 

5000 

5.275 ± .018 

5.275 ± .076 

1 

11000 

20000 

5.287 ± .006 

5.287 ± .024 

2 

18000 

80000 

5.291 ± .002 

5.291 ± .007 

3 

27000 

320000 
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Error reduction during mesh adaptivity Adapted mesh (18000 elements) 


Primal-Dual Problems in Fluid Mechanics 



Space- 

Time 

DG 

Tim Barth 

Introduction 

Cylinder 

Flow 

Nonlinear 

Conserva- 

tion 

Laws 

Space-time 

DG 

Error Rep- 
resentation 

Scalar 

transport 

Example 

Dual 

Problems 

Navier- 

Stokes 

Formulation 

Periodic 

Cylinder 

Flow 

Concluding 

Remarks 



Primal Mach number Dual x-momentum Adapted Mesh 



Software Implementation and extension to the Navier-Stokes 
Eqns 
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Space-Time FEM: 

o DG extension to the compressible Navier-Stokes equations using the 
symmetric interior penalty method of Douglas and Dupont, 1976) as 
described in Hartmann and Houston (2006) 

o Implements the discontinuous Galerkin discretization in entropy 
variables. 

o Unconditionally stable for all time step sizes 
o Solves both the primal numerical problem and the jacobian linearized 
dual problem arising in space-time error estimation, 
o High-order accuracy demonstrated in both ES and 


Dual Problems for Time Dependent Problems 
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Computing dual (backwards in time) problems looks expensive in terms of 
both storage and computation 

o Storage of the primal time slices for use in the locally linearized dual 
problem. 

o Approximation of the infinite-dimensional dual problem for the 
backwards in time dual problem. 

Tremendous simplification arising for periodic flow problems with period P 
when phase-independent functionals are utilized, e.g. mean drag 

o Functional independent of the startup transient 

o Only a small number of periods of the primal problem need be stored or 
recreated. 


□ i >00,0 


Periodic Cylinder Flow 
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Cylinder flow at Mach = 0.10, logarithm of |vorticity| contours 




Re=300 Re=1000 

Task: Represent and estimate the error in the mean drag force coefficient 

o Solve the primal problem using linear space-time elements 

o Construct a smooth phase invariant functional measuring the mean 
drag force coefficient 

o Solve the dual (backwards in time) problem using quadratic space-time 
elements 

o Calculated the estimated functional error and compare with a reference 
calculation using cubic elements 

< g ► <!► <!► 1 -0 0,0 


Mean Drag for Cylinder Flow 
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Dual problem, (f > ( x - rnom ) Dual defect, . 


Jdrag(u) = So Jr wall ( Force ■ tdrag) V(t)dx dt 

Example: Cylinder flow at Re=300 


< □ M 9 M | M | ► 1 -0 0,0 


Mean Drag Dual Problems at Re=300 and Re=1000 
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Dual problem at Re=300 
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Mean Drag for Cylinder Flow at Re=1000 
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Error representation buildup during the backward in time dual integration 


Error representation formula I 

Error indicator formula 




□ i -o o 


Mean Drag for Cylinder Flow at Re=1000 
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Adapted mesh from element indicators averaged over a period P 




Coarse mesh (1 2K elements) 2 level refined mesh (20K elements) 
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Non-Periodic Cylinder 
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Cylinder flow at Re=3900 and Re=10000 using quartic (p = 4) space-time 
elements. 

o Choosing measurement problems that are not genuinely stationary 
produces rapidly growing dual problems and dependency on the initial 
data. 



Dual solution corresponds to the average drag force over 3 approximate 
“periods”. 



Re=3900 



Re=1 000Gb 




Growth of Dual Problems 
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Reynolds Number 


Growth of drag functional dual solution 0 with increasing Reynolds number 


A Closing Note on the Use of Dual Problems in Uncertainty 
Quantification 


Space- 

Time 

DG 

Tim Barth 


Developing a capability to numerically compute primal/dual problems for 
compressible Navier-Stokes is a major undertaking. 
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Can this capability be reused in uncertainty quantification? 

Estep and Neckels (2006) observed that dual problems can be used to build 
a piecewise linear response surface for use in Monte Carlo (MC) and Quasi 
Monte Carlo (QMC) sampling of uncertain outputs when the output of 
interest is a functional. 
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J(u;u) : R m x R n — > R 

calculate statistics of the functional such as expectation 


/ J(u;cj) pdf ( uj) die = / J(u-,bj(fi))d n 
J v Jo 


V[J] = E[J 2 ] - E[J } 2 

using Monte Carlo (MC) or Quasi Monte Carlo (QMC) sampling. 


I -OH.O 


Fast MC and QMC Evaluation of Uncertain Output Functionals 
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Higher Order Parameter Sampling (HOPS), Estep and Neckels (2006) 

o Convert the statistics integration problem to uniform MC sampling on a unit 
hypercube. 

Partition the unit hypercube into smaller hypercube subdomains with size 
determined from accuracy of the linearized sampling formula. 

© In each hypercube subdomain C, center, calculate the primal solution u, and 
adjoint solution $,■ 

and the reduced sensitivity gradients (cf. A. Jameson, 1988) 
j dJ . . j dL 

9 = -7T {*,“/) 

OUJ OUJ 

© Apply MC or QMC integration in each C, using the linearized sampling formula 
for fixed values of J(u, , wf) and g J 

J(u, uj) m J(u„ wj) + gj (lj - u>i) 

< □ ► < g ► < 1 ► < 1 ► 1 -0 0,0 


Adaptive HOPS Surface 
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Estep and Neckels then consider adaptive refinement to improve approximation 
properties of the HOPS surface. 
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Original HOPS surface 



Adaptively refined HOPS surface 



< □ ► < if > 


>0 0.0 


Concluding Remarks 
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o Estimating and controlling numerical error in time-dependent 
calculations is fraught with difficulties 
o growth in backward-in-time dual problems, 
o loss of sharpness in error bounds. 

o The calculation of dual problems is computationally demanding 
o storage of primal time slices, 
o higher order solves of dual problem 

o Error representation/estimation results presented today barely scratch 
the surface 

o error control for general transient problems, 
o dual problems in the presence of flow bifurcations. 


< □ M 9 M | M | ► 1 -O Q, o 
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An Application of Error Estimation and Adaptive Error Control 
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Example: Euler flow past multi-element airfoil geometry. M = .1 , 5° AOA. 


lift coefficient 
(error representation) 

lift coefficient 
(error control) 

refinement 

level 

# elements 

equivalent uniform 
refinement 
# elements 

5.156 ± .147 

5.156 ± .346 

0 

5000 

5000 

5.275 ± .018 

5.275 ± .076 

1 

11000 

20000 

5.287 ± .006 

5.287 ± .024 

2 

18000 

80000 

5.291 ± .002 

5.291 ± .007 

3 

27000 

320000 



Error reduction during mesh adaptivity Adapted mesh (18000 elements) 



Example: Ringleb Flow 
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Schematic of Ringleb flow 
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Discontinuous Galerkin 


Iso-Density contours 


Example: A Scalar Time-Dependent PDE 
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Circular transport, A = (y, -x), of bump data 


ut + A ■ V u = 0 , 

u(x, 0) = V(1 / 1 0; x - x 0 ) , 


x e [-1,1]" 
x 0 = (7/10,0,0) 




Convergence, ||u - u h || /. 2 (n x [o,T]) 
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